Subroutine SIMULATION(EXP_L_prime,B, RET_PRIME, jExperimentNumber )

Use PARAMETERS
USE Functions_libs
USE Interpolate_func
USE random

Implicit NONE

Integer :: jtype1_orig, jht_new, jht_orig, jzt_int, jo_lagged, juntreated, joption_pay, jr_rep, jtype2_new, jdp_new, jdu_new, js_new,  I_pay,  jhc_ret, jmc, jtrpay, I_treat, jos, jtype4, jtype4_p, jno, jBenchmark_exp, jh_alt, jr_alt, jr_p_alt, jexp_new, jMEcat_p, jMEcat, jExperimentNumber , jspouse,counter, jzhc_p,jzhc, jmar,jmar_p, jtype3, jtype3_p, jzt, jzt_p,jsample,jas,jas_p,jsav,je,jh,jh_p, jr,jr_p,ji, jage,j, jht, jyears, jfix, jtype1,jtype2,jtype2_p,js,jdp,jdu,jx, jo,jo_p, jins, jhrs , jx_age, jo_act , jsav_act, kk  , jsav_upper , TR_indic, jinv_upper  , jsav_upper_ret  , jhk, jinc

Real (Kind=dd), dimension (n_as,n_type1,n_type3,n_hc_grid,n_ret-1,n_o,n_mar,2)::  EXP_L_prime  !value functions    
Real (Kind=dd), dimension (n_o,n_type4,n_zt,n_as,n_type3,n_type1,n_hc_grid,n_ret-1)::  B !value function while working. (1=not have option to treat and not pay, 2=have option)

Real (Kind=dd), dimension (n_h,n_r,n_o,n_mar,n_zhc)::  EXP_L_prime_INT

Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_mar):: RET !value functions when retired
Real (Kind=dd), dimension (n_hc_ret,n_MEcat,n_as,n_type1,n_h,n_r,n_type2,n_ret : n_age,n_mar):: RET_PRIME !value functions when retired
Real (Kind=dd), dimension (n_MEcat,n_h,n_r,n_type2,n_mar):: RET_PRIME_INT

Real  (Kind=dd) :: Resources, prob_s, prob_d, income_bt, income, base_1, base_2 , Total_tax, actual_earn, tax_SS, tax_med, V_death, V_c_all_trpay, V_cf_trpay
Real (Kind=dd) :: last_value, maxB, valret   , check
Real (Kind=dd) :: UTI,x, EXPECTED_RET_PRIME 

Real (Kind=dd), dimension (1:40)::  zt_hypo
Real (Kind=dd), dimension (1:40):: o_act_bench, sav_bench, HC_bench, wage_bench, o_offer_hypo

Real (Kind=dd), dimension (1:443)::  ran_hypo ! these are the saved random draws from the simulation without shocks. - fill in those with shocks. 6+11*39+5 (6 extra for age 1, and at age 40, we only need to draw 5 shocks that are relevant for decisions at age 40 (64). All other shocks are not relevant until age 65).

Real (Kind=dd):: CON, TR, ASP , MAX_SAV  , as_opt, TR_opt ,   INT_sav_upper_ret   , INT_sav_upper, EXPECTED_B_PRIME
Real (Kind=dd), dimension (n_o) :: EXPECTED_B

Real (Kind=dd) :: jx_int_p, prob_sim
Real (Kind=dd), dimension (n_zhc) :: HK_INTERPOL_sim

Real (Kind=dd):: Z_INT_ALT , Z_INT_OFFER_ALT , temp,  jcons_act, jas_p_act, jtr_act  ,jassets  , jx_int, Z_INT , Z_INT_OFFER , B_act , RET_act  , jx_age_interpol  , EXPECTED_RET_PRIME_act, EXPECTED_B_PRIME_act
Integer :: bb,d,o,f   , gr, gg
Real (Kind=dd) :: rmaxipzt  , rmaxipzhc
Real (Kind=dd), dimension (2,n_MEcat,n_type2,n_o) ::  B_int

Real (Kind=dd), dimension(3)::  CON_sim,  ASP_sim, B_sim,   TR_sim
Real (Kind=dd), dimension(2):: income_pay, Total_tax_pay, base_2_pay, base_1_pay
                    
Integer :: indicator_ret, idum, ind_death, intepsi
Real (Kind=dd), dimension(n_mar) :: rmaxiTPzp
Real (Kind=dd) ::  rmaxiIDzp
Real (Kind=dd) :: rx, prob, ry, rb, rd, rf, ri , ro, rinv, i
Real (Kind=dd) :: epsi,weight_epsi, zp_step

1008 Format(I5,1x,I2, 1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,I1,1x,F6.2 , 1x, F6.2 , 1x, F9.2 , 1x,F9.2, 1x,F9.2, 1x,F10.2, 1x,F10.2, 1x, F9.2, 1x, F8.2 , 1x, F10.2, 1x, F4.1, 1x,F9.2, 1x, F10.2,1x,I1,1x, F9.2,1x,F10.2,1x,F10.2,1x,F10.2,1x,F10.2, 1x, I1, 1x, I1, 1x, I1, 1x, I1, 1x, I1 , 1x, F10.6)
      ! file for writing hk - write for all ages from 30 onwards, 40, 50, 60 if no shocks at these ages and no deterioration in health
1004  Format(F16.11)  ! this is for human capital
1005 Format(F16.11) ! this is for random shocks 
1006 Format(I1)  ! zt shocks are 1, 2, or 3
     
     
!Grid for savings is finer in the simulation     
sav_sim(:)=0.0d0
open(258, file='External_Parameters/savings_grid_sim.txt')
do jsav=1, n_sav_sim
   read (258,*) sav_sim(jsav)
end do 
close (258)
sav_sim(:)=sav_sim(:)/5200.0d0     

! sav_grid: useful for interpolation to find the maximum savings
sav_grid_sim(:)=0.0d0
 do jsav=1, n_sav_sim
 sav_grid_sim(jsav) = dble(jsav)
 end do
     
! NUMBER OF PEOPLE OF EACH TYPE SIMULATED
! I simulate histories of n_sample agents. In_type1(jtype1)  tells how many agents of each type to be simulated. (There are 18 types in total). 
do jtype1=1, n_type1
n_samp(jtype1)= dble(n_sample)*In_type1(jtype1)  
end do







!***********************************************************************
! Simulating the economy
!***********************************************************************

! jx_age is an integer number on a constructed grid for human capital and age, e.g., 1 (age=1, hk=1), 2 (age=2, hk=1), 3 (age=2, hk=2) etc
! jx_age_interpol is a real number that can take values in between those of jx_age.  E.g., 4.5 (age=3, hk=1.5, meaning .5 yrs of full time work)
  
 do jBenchmark_exp= 21, 24 ! we also need to do 28
 
 if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0 ) then
open(12, file='Simulated_Data/Simulated_data0.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')
else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.5 ) then
open(12, file='Simulated_Data/Simulated_data05.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21 ) then
open(12, file='Simulated_Data/Simulated_data021.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')
open(16, file='Simulated_Data/Simulated_wages.txt')
open(17, file='Simulated_Data/Simulated_o_act.txt')
open(18, file='Simulated_Data/Simulated_hum_cap.txt')
open(19, file='Simulated_Data/Simulated_H.txt')
open(20, file='Simulated_Data/Simulated_savings.txt')
open(21, file='Simulated_Data/Simulated_zt.txt')
open(22, file='Simulated_Data/Random_numbers_old0.txt')
open(23, file='Simulated_Data/Simulated_o_offers.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.22 ) then ! same o_act and savings as benchmark
open(12, file='Simulated_Data/Simulated_data022.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')
open(17, file='Simulated_Data/Simulated_o_act.txt')
open(20, file='Simulated_Data/Simulated_savings.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.23 ) then
open(12, file='Simulated_Data/Simulated_data023.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')
open(18, file='Simulated_Data/Simulated_hum_cap.txt')
open(21, file='Simulated_Data/Simulated_zt.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.24 ) then
open(12, file='Simulated_Data/Simulated_data024.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')
open(16, file='Simulated_Data/Simulated_wages.txt')

else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.28 ) then
open(12, file='Simulated_Data/Simulated_data028.txt')
open(14, file='Simulated_Data/Random_numbers0.txt')
open(22, file='Simulated_Data/Random_numbers_old0.txt')
open(23, file='Simulated_Data/Simulated_o_offers.txt')
open(18, file='Simulated_Data/Simulated_hum_cap.txt')
open(21, file='Simulated_Data/Simulated_zt.txt')

end if 

do jfix=1, min(n_fix, I_state_fix) 
    je = 1 
    do jht=1, min(n_ht, I_state_ht)
    jht_orig=jht
    jht_new=jht
    
    call TYPE1_DET(je,jht_new,jfix,jtype1)
    call TYPE1_DET(je,jht_orig,jfix,jtype1_orig)
        
do jsample=1,int(n_samp(jtype1_orig))   
    
    ! read the same random draws as benchmark
      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.ge.21) then
        do jno=1, 443
            read (14,*) ran_hypo(jno) ! 
        end do
      end if 
      
    if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.22) then ! read LS and savings
        do jno=1, 40
            read (17,*) o_act_bench(jno)
            read (20,*) sav_bench(jno)
        end do    
     else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.23) then ! read HC
        do jno=1, 40
            read (18,*) HC_bench(jno)
            read (21,*) zt_hypo(jno)
        end do
    else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.24) then ! read wages
        do jno=1, 40
            read (16,*) wage_bench(jno)
        end do
    else if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.28) then ! read HC
        do jno=1, 40
            read (18,*) HC_bench(jno)
            read (21,*) zt_hypo(jno)
            read (23,*) o_offer_hypo(jno)  
        end do
    end if 
    
do jage=1,n_age 

! *************   FIRST PERIOD WHEN AGE=1    **************
if (jage.eq.1) then 
               jassets=assets_initial(je)
               jx=0 ! 0 years of experience
               jx_int=0.d0    

!Pick jh from the initial distribution of health at age 25 from data
                    rb=ran0(idum) 
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                        write(14,1005) rb 
                        else 
                             rb=ran_hypo(1) 
                         end if 
                    bb=1
                    prob=In_h(je,bb,jht_orig)
                 57 if (rb.lt.prob) then
                    jh=bb
                    else 
                    bb=bb+1
                    if (bb.le.n_h) then
		            prob=prob+In_h(je,bb,jht_orig)
		            go to 57
		            else 
		            jh=n_h ! I added this because due to rounding up errors
                    end if 
                    end if    
                    
                      jh_alt=jh
                    
!Pick jr from the initial distribution of health at age 25 from data
                    rb=ran0(idum) 
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb
                          else
                           rb=ran_hypo(2) 
                       end if 
                    bb=1
                    prob=In_r(je,bb)
                 58 if (rb.lt.prob) then
                    jr=bb
                    else 
                    bb=bb+1
                       if (bb.le.n_r) then
		            prob=prob+In_r(je,bb)
		            go to 58
                    else 
		            jr=n_r ! I added this because due to rounding up errors
                    end if 
                    end if 
                    
                    jr_alt=jr 
                    
!Pick jmar from the initial distribution at age 25 from data
                    rb=ran0(idum)   
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(3) 
                       end if 
                    prob=In_mar(je,jh)
                    if (rb.le.prob) then
                    jmar=2
                    else 
                    jmar=1
                    end if
                    
!draw spouse work status
                        rb=ran0(idum) 
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(4) 
                       end if 
                    
                    if (jmar.eq.1) then
                        jos=1 ! assign any number, does not matter
                    else if (jmar.eq.2) then ! draw spouse work status
                        if (rb.lt.(prob_os(jage,je,jh,1,jfix))) then
                        jos=1
                        else 
                        jos=2
                        end if 
                    end if 
                    
!Draw employment offer
                  rb=ran0(idum) ! draw a random # between 0-1
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(5) 
                          end if 
                        
                           if  (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.28) then       
                           jo= o_offer_hypo(jage)
                           else 
		            bb=1
		            prob=Pr_o_25(bb, je)
85		            if (rb.lt.prob) then
		            jo=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_o) then
		            prob=prob+Pr_o_25(bb, je)
		            go to 85
                    else 
                        jo=n_o
		            end if
                    end if
                    end if
		            
!Draw zt 
                    rb=ran0(idum) ! draw a random # between 0-1
                    
                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                          write(14,1005) rb 
                          else
                           rb=ran_hypo(6) 
                          end if 
                     
                 if  (jExperimentNumber.EQ.0.AND.((jBenchmark_exp.EQ.23).OR.(jBenchmark_exp.EQ.28))) then       
                           jzt= zt_hypo(jage)
                           else 
		            bb=1
		            prob=p_zt(je,bb,2)
86		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,2)
		            go to 86
                    else 
                        jzt=n_zt
		            end if
                    end if 
                 end if 
                    
end if  !for jage=1         
!*********************************************************

! *************   ALL AGES<N_RET    **************
if (jage.lt.n_ret) then	
    
    !experiments 
        if (jExperimentNumber.EQ.0.AND.((jBenchmark_exp.EQ.23).OR.(jBenchmark_exp.EQ.28)).AND.(HC_bench(jage).ne.100.0d0)) then
        jx_int = HC_bench(jage)
        jzt= zt_hypo(jage)
        end if 
        
        if (jExperimentNumber.EQ.0.AND.(jBenchmark_exp.EQ.28).AND.(HC_bench(jage).ne.100.0d0)) then
        jo=o_offer_hypo(jage)
        end if 
    
    !Draw mortality shock            
                     rb=ran0(idum)
                                 if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                                write(14,1005) rb 
                                end if 
                                 if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.ge.21)  then
                                rb=ran_hypo(7+11*(jage-1))
                                    end if
		             if (rb.gt.p_surv(jh,je,jmar,  jage)) then !individual dies	 
     !filling in random draws to files                     
if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
    if (jage.eq.40) then
        do jno=1,7
            rb=ran0(idum)
            write(14,1005) rb
        end do
        
        else 
                 do jno=1, 18 +(39-jage)*11 ! this makes sure that all individuals have 443 observations written to the random draws file- makes it easier to read.  
                     rb=ran0(idum)
            write(14,1005) rb
                 end do 
        end if 
end if  

if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then ! write implausible numbers so we know
 
                  do jno=1,(41-jage)   
                        write (17,1006) 6 !o_act_bench
                        write (20, 1005) 100.0d0 !sav_bench 
                        write (18,1005)  100.0d0 !HC_bench 
                        write (16,1005) 0.0d0 ! wage_bench ! wage offer
                        write (19,1006) 0 !H
                        write (21,1006) 0 !zt_hypo
                        write (23,1006) 6 !o offer benchmark
                  end do 
                  
               ! now write old people shocks   
               do jno=1, 213   
                     rb=ran0(idum)
            write(22,1005) rb
                 end do 
      
end if

                     go to 33 !simulate new individual 
                     end if

call TYPE4_DET(jmar,jos,jtype4) ! if jage is not 1, then jmar and jos have been already drawn towards the end of the code.                                            
call TYPE3_DET(jh,jr,jtype3) ! we need this index only for non-retirees

! determine optimal employment status "jo_act"
   !The individual decides whether to accept the offer or decline
   !interpolate with respect to jas and jx.

if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.22.AND.(o_act_bench(jage).ne.0)) then !the last condition identifies cases where people had died in the benchmark
    jo_act=o_act_bench(jage)
    else 

 if (jo.eq.1) then 
         jo_act=1  ! THIS IS YOUR ACTUAL EMPLOYMENT STATUS 
         
 else  ! HAVE TO DECIDE BETWEEN WORKING AND NOT WORKING
         EXPECTED_B(jo)=0.d0
         EXPECTED_B(1)=0.0d0

call interpol2d(as(1:n_as,jtype1,jage), xfn(   1: n_hc_grid , jage)  ,   B(1,jtype4,jzt,1:n_as,jtype3,jtype1,1: n_hc_grid,jage) , jassets, jx_int, EXPECTED_B(1))   

call interpol2d(as(1:n_as,jtype1,jage), xfn(   1: n_hc_grid , jage)  ,  B(jo,jtype4,jzt,1:n_as,jtype3,jtype1,1: n_hc_grid,jage) , jassets, jx_int, EXPECTED_B(jo))              
         If   (EXPECTED_B(jo).gt.EXPECTED_B(1)) then  
           jo_act=jo   
         else 
           jo_act=1 
         end if 
                
 end if 
 end if 
!Draw health shocks (TYPE2)
                  rb=ran0(idum) ! draw a random # between 0-1 
                  
            if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
            write(14,1005) rb
            else 
                rb=ran_hypo(8+11*(jage-1))
            end if
                 
		            bb=1
		            prob=shocks_risk(je,bb, jh, jr, jage) 
		         90 if (rb.lt.prob) then
		            jtype2=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_type2) then
		            prob=prob+shocks_risk(je,bb, jh, jr, jage)
		            go to 90
		            else 
		            jtype2=n_type2 !go to 89  ! I added this because due to rounding up errors
		            end if
                    end if

if ((jExperimentNumber.EQ.0).AND.((jBenchmark_exp.EQ.5).OR.(jBenchmark_exp.EQ.22).OR.(jBenchmark_exp.eq.23).OR.(jBenchmark_exp.eq.24).OR.(jBenchmark_exp.eq.28) )) then
 jtype2=1
end if	                    
                    
if (jtype2.eq.1) then
jdu=1 ; jdp=1 ; js=1
else if  (jtype2.eq.2) then
jdu=1 ; jdp=2 ; js=1
else if  (jtype2.eq.3) then
jdu=1 ; jdp=1 ; js=2
else if  (jtype2.eq.4) then
jdu=1 ; jdp=2 ; js=2
else if  (jtype2.eq.5) then
jdu=2 ; jdp=1 ; js=1
else if  (jtype2.eq.6) then
jdu=2 ; jdp=2 ; js=1
else if  (jtype2.eq.7) then
jdu=2 ; jdp=1 ; js=2
else if  (jtype2.eq.8) then
jdu=2 ; jdp=2 ; js=2
end if
		             
!draw the medical expenditure catastrophic shock
                    rb=ran0(idum) ! draw a random # between 0-1
                    
                                if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
            write(14,1005) rb 
            else
        rb=ran_hypo(9+11*(jage-1))
                      end if
                                     
		            bb=1
		            prob= Prob_MEcat(bb)
		         28 if (rb.lt.prob) then
		            jMEcat=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_MEcat) then
		            prob=prob+Prob_MEcat(bb)
		            go to 28
		            else 
		            jMEcat= n_MEcat  !go to 27  ! I added this because due to rounding up errors
		            end if
                    end if  
                    
                    
                    ! draw whether you have the option to not pay and treat. 
                                        rb=ran0(idum) ! draw a random # between 0-1
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
            write(14,1005) rb 
            else
        rb=ran_hypo(10+11*(jage-1))
                      end if
                                        
       
                    
            		prob= Shock_type(1,jo_act,jh,jage) + Shock_type(2,jo_act,jh,jage)   
                  if (rb.le.prob) then
		            joption_pay=1 ! cannot treat and not pay
		            else 
		            joption_pay=2 ! can treat and not pay
                    end if  
                  
                    
                    
                    ! if joption_pay=1, then if you are uninsured, with some probability you will not treat.
                                        rb=ran0(idum)
                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0) then
                                            write(14,1005) rb 
                                            else
                         rb=ran_hypo(11+11*(jage-1))
                      end if
                      
                    if ((joption_pay.eq.1).and.(jtype2.ne.1)) then
                    if (rb.lt.(Shock_type(1,jo_act,jh,jage)/(Shock_type(1,jo_act,jh,jage)+Shock_type(2,jo_act,jh,jage)))) then ! you cannot treat 
                        juntreated=1 ! cannot treat
                    else 
                        juntreated=2 ! this means you have option to treat
                    end if 
                    else 
                        juntreated=2 ! set to 2 for everyone else 
                    end if 
                    
      ! we only care about the charges for uninsured. so we actually use the OOP bc for uninsured this is MC*.6              
         if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(1)) then 
        jmc=1
        else if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(2)) then
        jmc=2
        else if (oop(jtype2,jh,jage,jMEcat,jo_act).le.MC_TS(3)) then
        jmc=3
        else
        jmc=4
        end if
 
!given the optimal work decision and sick days, determine human capital at the beginning of next period, BEFORE the human capital shock occurs. "jx_int_p"             

        if ((jExperimentNumber.EQ.0).AND.((jBenchmark_exp.EQ.23).OR.(jBenchmark_exp.EQ.28)).AND.(HC_bench(jage+1).NE.100.0d0)) then
                                jx_int_p= HC_bench(jage+1)
                                   else
                                    
    if (jo_act.eq.1) then 
     jx_int_p=jx_int
     else if (jo_act.eq.2.OR.jo_act.eq.3) then  
     jx_int_p= min((jx_int+ (max( (hrs_pt - phi(je,jh,jtype2) ), 0.0d0) /hrs_ft )),  dble(n_hc) )
     else if (jo_act.eq.4.OR.jo_act.eq.5) then 
     jx_int_p=   min((jx_int+ (max((hrs_ft - phi(je,jh,jtype2) ), 0.0d0) /hrs_ft )),  dble(n_hc) ) ! jx_int goes up by 1 if you work full time and have no sick days. We cap it at 100. 
     end if
     
     end if 
     
!consider each possible shock to human capital and get the human capital levels next period associated with each. - note probabilities do not come into play yet. "HK_INTERPOL_sim(jzhc_p)"   
if (jage.le.n_ret-2) then    
do jzhc_p=1,n_zhc ! allow 3 possible shocks
      HK_INTERPOL_sim(jzhc_p) =max(0.0d0,  min (jx_int_p *zhc(jzhc_p), dble(n_hc))) !restrict between 0-100
end do       
end if

!Find wages associated with current human capital level "jx_int"-- NEED TO INTERPOLATE: get "Z_INT" and "Z_INT_OFFER"
Z_INT_OFFER=0.0d0
Z_INT=0.0d0

 if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.24.AND.(wage_bench(jage).NE.0.0d0)) then
        Z_INT_OFFER=wage_bench(jage)
        else

call InterpolateScalar(HK(0 : n_hc), z(0 : n_hc,jh,je,jfix,jzt), jx_int, Z_INT_OFFER)
        end if 
        
    if (jo_act.ge.2) then  
     Z_INT=Z_INT_OFFER
    end if
    
if (jage.eq.(n_ret-1)) then 
!Determine human capital group for the purpose of SS next period for those aged 64 (jhc_ret)
if (jx_int.le.hc_TS(1,je)) then 
    jhc_ret=1
else if (jx_int.le.hc_TS(2,je)) then 
    jhc_ret=2
else 
    jhc_ret=3
end if
else 
   jhc_ret=1  ! just to have a value although this will not matter for anything. 
end if 
 

! calculate incomes if pay or not pay
 actual_earn=0.0d0
 income_bt=0.0d0
 tax_SS=0.0d0
 tax_med=0.0d0
 
  actual_earn = max(0.0d0 , ( wage_penalty(jo_act,je)*Z_INT * (hrs_offer(jo_act) -phi(je,jh,jtype2) )) ) ! earnings net of sick days
  income_bt= ((interest-1.0d0))*jassets +  actual_earn  ! income before tax is equal to interest income plus earnings net of sick days
 
if ((jo_act.eq.3).OR.(jo_act.eq.5)) then ! husband holds insurance -- so deduct from his income only
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn -prem_ephi(jo_act,jtype4) ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix)  )) !! own plus spouse's
tax_med = tau_med* max(0.0d0,  ( actual_earn -prem_ephi(jo_act,jtype4))) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) ))
else ! subtract from wife. note this can be zero
tax_SS = tau_ss*max(0.0d0, min(yhat ,  actual_earn  ))   + tau_ss*max(0.0d0, min(yhat ,  inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) )) !! own plus spouse's
tax_med = tau_med* max(0.0d0,  ( actual_earn )) +  tau_med* max(0.0d0,  ( inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)))
end if   
  
 ! for all of these, pay OOP =1, not pay OOP=2 
 income_pay(:)=0.0d0
 Total_tax_pay(:)=0.0d0
 base_2_pay(:)=0.0d0
 base_1_pay(:)=0.0d0
 
        ! if pay OOP 
  income_pay(1) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4) -tax_SS - tax_med - max(0.d0, oop(jtype2,jh,jage,jMEcat,jo_act) + oop_spouse(jage,je,jtype4,jo_act) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(1) = max(0.0d0,  (( (income_pay(1)*5200.0d0) - tax_lambda(jmar) *( (income_pay(1)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(1) =  prem_ephi(jo_act,jtype4) +  oop(jtype2,jh,jage,jMEcat,jo_act) +  oop_spouse(jage,je,jtype4,jo_act) +Total_tax_pay(1)  ! total resources subtracted 
  base_1_pay(1) =  interest*jassets + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(1)        ! cash on hand (if you pay the OOP)
      ! if not pay OOP
  income_pay(2) = max(0.d0,  income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -prem_ephi(jo_act,jtype4)  -tax_SS - tax_med -  max(0.d0,  oop_spouse(jage,je,jtype4,jo_act) - 0.075d0*(income_bt+  inc_spouse(jage,je,jh,jtype4,jfix)) ) )  ! taxable income
  Total_tax_pay(2) = max(0.0d0,  (( (income_pay(2)*5200.0d0) - tax_lambda(jmar) *( (income_pay(2)*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) +  tax_SS + tax_med  ! tax you have to pay
  base_2_pay(2) =  prem_ephi(jo_act,jtype4) +  oop_spouse(jage,je,jtype4,jo_act) + Total_tax_pay(2)  ! total resources subtracted 
  base_1_pay(2) =  interest*jassets + inc_spouse(jage,je,jh,jtype4,jfix)  + actual_earn - base_2_pay(2)        ! cash on hand (if you do not pay the OOP)
  
!Determine income group
          if (income_bt.le.INC_TS(1)) then 
        jinc=1
        else if (income_bt.le.INC_TS(2)) then 
        jinc=2
        else if (income_bt.le.INC_TS(3)) then 
        jinc=3
        else if (income_bt.le.INC_TS(4)) then 
        jinc=4
        else 
        jinc=5
        end if
    
        
!DETERMINE IF THE PERSON WILL PAY/TREAT and OPTIMAL CONSUMPTION/SAV
!(determine assets next period): jcons_act, jas_p_act , jsav_act

  !initiate        
CON_sim(:)=0.0d0
  ASP_sim(:)=0.d0
  B_sim(:)=0.d0
  TR_sim(:)=0.0d0
  
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
  V_cf_trpay=0.0d0
  V_c_all_trpay=0.0d0
  V_death=0.0d0
        jcons_act=0.0d0
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
    if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.22.AND.(o_act_bench(jage).NE.0)) then
    jtrpay=1
    I_treat=2
    I_pay=1
    jas_p_act = sav_bench(jage)
    jcons_act = max( (income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1) - jas_p_act) , cons_floor_at(jage, jmar,je,jh))
    B_sim(1)=-100.0d0 ! just a number to fill in; don't need it
        if ( jcons_act .le. (cons_floor_at(jage, jmar,je,jh))) then 
          jtr_act=  cons_floor_at(jage, jmar,je,jh) - (income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1) - jas_p_act)  
        end if 
        
    else 
  
!If I don't have a shock, I pay medical bill. 
if (jtype2.eq.1) then
     jtrpay=1
    I_treat=2
    I_pay=1
    
    !If on cons floor
    if ( base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) then 
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(1)
    
    ! still need to calculate V because we need it all the time for welfare analysis
    ! get utility
    call func_uti(jage,1,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 1
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(1) = UTI + beta_hat(je)* EXPECTED_B_PRIME 

    else
              jtr_act=0.d0 
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(1)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if
               
!VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,1,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets   
                    ! the value function tomorrow depends on whether I'm working today
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,1,je,jh,jo,jtype4,jtype2,(base_1_pay(1)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(1) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,1,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
               else 				
					call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
               end if 
               
                  EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 

     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
51 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
        B_sim(1)=   maxB    
    end if 
    
 !If I have a shock   
else 
if (( base_1_pay(2) .le. (cons_floor_at(jage, jmar,je,jh))).AND.(juntreated.EQ.1).AND.(I_treat_Mcaid.EQ.0)) then ! if I'm on the cons floor even if I don't pay AND cannot treat

        jtrpay=3 ! you don't treat
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(2)
  
    ! get utility
    call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 3
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,1,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(2,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(jtrpay) = UTI + beta_hat(je)* EXPECTED_B_PRIME 
    
else if (( base_1_pay(2) .le. (cons_floor_at(jage, jmar,je,jh))).AND.((juntreated.EQ.2).OR.(I_treat_Mcaid.EQ.1))) then ! if I'm on the cons floor even if I don't pay AND can treat - then you always "treat and pay"
    jtrpay=1
    I_treat=2
    I_pay=1
    jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
    jas_p_act=0.d0 
    jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(1)
    
    ! get utility
    call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,jcons_act,UTI) ! second argument is trpay which we set to 1
             V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.           
     if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
     end if 
     
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
                     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         B_sim(jtrpay) = UTI + beta_hat(je)* EXPECTED_B_PRIME 
    
else
    
if (joption_pay.eq.2) then   ! if I have all 3 options 
do jtrpay=1, 3
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    
    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else if (jtrpay.eq.2) then
        I_treat=2 ! treat
        I_pay=2 ! not pay
    else 
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.     
                    if (jo_act.eq.1) then 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      else               
          EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      end if 
      
     end do ; end do ; end do   ; end do; end do  
  
         end if 
         
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets  
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                    end if 
                    
      if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat              
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
      else
                              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

      end if 
      
      end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 
 
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.   
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     
     else

EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                                        if (jo_act.eq.1) then

                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                                        else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                                        end if 
                                        
               if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat     
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
               else 
                                       EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

               end if 
               
               end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 

     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
end do ! jtrpay

! determine best option and reset I_treat and I_pay
if ((B_sim(1).eq.0.0d0).OR.(B_sim(2).eq.0.0d0).or.(B_sim(3).eq.0.0d0)) then
PRINT *,'Error'
pause
end if 
if ((B_sim(1).gt.B_sim(2)).AND.(B_sim(1).ge.B_sim(3))) then
    jtrpay = 1 
else if ((B_sim(2).ge.B_sim(1)).AND.(B_sim(2).ge.B_sim(3))) then
    jtrpay = 2 
else if ((B_sim(2).eq.B_sim(1)).AND.(B_sim(2).eq.B_sim(3))) then
    jtrpay = 2 
    else
    jtrpay = 3 
end if 

jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay)

 else if (((joption_pay.eq.1).AND.(juntreated.EQ.2)).OR.(((joption_pay.eq.1).AND.(juntreated.EQ.1)).AND.(I_treat_Mcaid.EQ.1).AND.(base_1_pay(1) .le. (cons_floor_at(jage, jmar,je,jh))) )) then ! if you don't have the option to treat and not pay. BUT HAVE OPTION TO TREAT

    do jtrpay=1, 3, 2 ! i only consider 1 and 3       
        
        ! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    
    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else 
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        !jsav_act=0.0d0 
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 
 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.  
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    
            
               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets    
                     if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                     else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                                             end if 
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
                             EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.  
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    end if 
                    
                        
          if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
          else 
                   EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

          end if 
          
          
          end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    else 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

                    end if 
                    
                    
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
    end do ! jtrpay
    
    
    ! determine best option and reset I_treat and I_pay
if ((B_sim(1).eq.0.0d0).or.(B_sim(3).eq.0.0d0)) then
PRINT *,'Error'
pause
end if 
if (B_sim(1).ge.B_sim(3)) then
    jtrpay = 1 
    else
    jtrpay = 3 
end if 

jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay)


 else if ((joption_pay.eq.1).AND.(juntreated.EQ.1)) then ! you cannot treat

    
    jtrpay=3 ! you don't treat
    
! I_treat and I_pay    
    
!jtrpay=1 	pay and get treated
!jtrpay=2	don't pay, get treated
!jtrpay=3	don't pay, not treated
    

        I_treat=1 ! not treat
        I_pay=2 ! not pay


!initiate
  CON=0.0d0
  ASP=0.d0
  UTI=0.d0
  TR=0.0d0
  MAX_SAV=0.0d0 
        jcons_act=0.0d0
        !jsav_act=0.0d0 
        jas_p_act=0.0d0
        jtr_act=0.0d0
        B_act=0.0d0
        
!If you are on the consumption floor        
if ( base_1_pay(I_pay) .le. (cons_floor_at(jage, jmar,je,jh))) then ! you get transfer. you only get to this point if you are not on transfer if you don't pay. so you'll still have to calculate the value if you don't pay
 jcons_act=cons_floor(jage, jmar,je,jh)  ! cons floor is at the HH level; CON is after tax HH consumption - it's what goes into utility fn
 jas_p_act=0.d0 
 jtr_act=  cons_floor_at(jage, jmar,je,jh) - base_1_pay(I_pay) 
 call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI)
 V_death =  Bequest_val(0.0d0)  + cost_death 
 
!continuation value
 EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then 


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.   
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
     if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
              EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
maxB=UTI+ beta_hat(je)*EXPECTED_B_PRIME 
 
else
  jtr_act=0.d0
! find maximum savings on the grid (jsav_upper)
  MAX_SAV= income_bt + inc_spouse(jage,je,jh,jtype4,jfix) -base_2_pay(I_pay)  - cons_floor_at(jage, jmar,je,jh)   
  
  !find lower bound 
          jsav_upper=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper+1) )
           jsav_upper=jsav_upper +1
           if (jsav_upper.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper.lt.1) then
                 PRINT *,'Error'
               end if 
                    

               
  !VALUE ASSOCIATED WITH CONSUMING on the consumption floor and saving max amount: V_cf_trpay
         call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,cons_floor(jage, jmar,je,jh),UTI) 
         V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    RET_PRIME_INT(:,:,:,:,:)=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                    !interpolate for assets
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), (MAX_SAV+jassets), RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )  
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital and assets 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), (MAX_SAV+jassets), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    end if 
                    
                    
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
     else 
                             EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

     end if 
     
     
     end do ; end do ; end do   ; end do; end do  
  
         end if 
                 
         V_cf_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME  

         
 !VALUE ASSOCIATED WITH CONSUMING ALL AND 0 ASSETS: V_c_all_trpay
         
           call func_uti(jage,jtrpay,je,jh,jo,jtype4,jtype2,(base_1_pay(I_pay)/(1.0d0+taxc)),UTI) 
         V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         EXPECTED_B_PRIME=0.0d0
         
         if (jage.eq.(n_ret-1)) then 

                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
      !interpolate with respect to human capital. Assets will be 0 next period.     
                    if (jo_act.eq.1) then
     call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )
                    else 
                             call InterpolateScalar(xfn(  1: n_hc_grid , jage+1) , EXP_L_prime(1,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) )

                    end if 
                    
          if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
     EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))
          else 
                   EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (p_surv(jh_p,je,jmar_p, jage+1)*EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p) + V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)))

          end if 
          
          
          end do ; end do ; end do   ; end do; end do  
  
          end if 
                 
         V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
         
          
!Iterate over possible savings and find optimal   
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
do jsav= jsav_upper , 1, -1
! get ASP, CON, and UTI associated with this jsav    
        ASP=  jassets+  sav_sim(jsav)
     if (ASP.le.0.d0)  exit   ! I already did the case where you consume all and save 0, so that's why we have "le," not "lt."
        V_death= cost_death + Bequest_val(ASP)  
        CON= ( income_bt + inc_spouse(jage,je,jh,jtype4,jfix) - base_2_pay(I_pay) -sav_sim(jsav) )/(1.0d0+taxc)
        call  func_uti(jage,jtrpay,je,jh,jo_act,jtype4,jtype2,CON,UTI)           
 
!Solve continuation value given this jsav and jtrpay
          EXPECTED_B_PRIME=0.0d0
         if (jage.eq.(n_ret-1)) then !need to interpolate only for assets; for human capital, we discretize it into bigger categories, so no need to interpolate.
                    RET_PRIME_INT(:,:,:,:,:)=0.0d0


                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                        call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME +Tp_h(jh_p, jh,je,jht_new,jage, jtype2,I_treat,jmc)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p)* p_surv(jh_p,je,jmar_p, jage+1) +  V_death * (1.0-p_surv(jh_p, je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do
  
         else 
             EXP_L_prime_INT(:,:, :,:,:)=0.0d0
                    do jo_p=1, n_o  ; do jh_p=1, n_h ; do jr_p=1, n_r ;  do jmar_p=1, n_mar;  do jzhc_p=1,n_zhc
                    call TYPE3_DET(jh_p,jr_p,jtype3_p) 
                    if (jo_act.eq.1) then
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,1), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           
                    else 
                    call interpol2d(as(1:n_as,jtype1,jage+1), xfn(   1: n_hc_grid , jage+1), EXP_L_prime(1:n_as,jtype1,jtype3_p,1: n_hc_grid,jage+1,jo_p,jmar_p,2), ASP, HK_INTERPOL_sim(jzhc_p), EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p))                           

                    end if 
                    
                                            
                    if ((jtrpay.eq.1).OR.(jtrpay.eq.2)) then ! treat 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 
                    else 
                    EXPECTED_B_PRIME= EXPECTED_B_PRIME +  prob_s_treat_new(jmc+1,jage, jo_act,je,jht_new,jtype3,jtype2,jo_p,jtype3_p,jzhc_p)*Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act)  *(EXP_L_prime_INT(jh_p,jr_p, jo_p,jmar_p,jzhc_p)* p_surv(jh_p,je,jmar_p, jage+1) + V_death * (1.0-p_surv(jh_p,je,jmar_p, jage+1)) ) 

                    end if 
                    
                    
                    end do ; end do ; end do   ; end do; end do  
  
          end if 
                                 
B_act=UTI+ beta_hat(je)*EXPECTED_B_PRIME 


     if (B_act.lt.maxB) exit
              maxB=B_act  
              jas_p_act=ASP
              jcons_act=CON
    
       
 end do !jsav
! Found jcons_act , jas_p_act and B_act assocaited with optimal savings. 
   
! now also need to compare with:    V_c_all_trpay. And if this is max, then set jas and jcons accordingly.
  if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= base_1_pay(I_pay)/(1.0d0+taxc)
  end if 
  

end if ! for being on the consumption floor

! save CONS, ASP, TR and B for the 3 jtrpay options. 

  CON_sim(jtrpay)=jcons_act
  ASP_sim(jtrpay)=jas_p_act
  TR_sim(jtrpay)=jtr_act
  B_sim(jtrpay)=maxB
  
  
  
jtr_act= TR_sim(jtrpay)
jas_p_act=ASP_sim(jtrpay)
jcons_act=CON_sim(jtrpay) 
    
    

    
    
end if ! for having the options in trpay


    if (jtrpay.eq.1) then
        I_treat=2 ! treat
        I_pay=1 ! pay
    else if (jtrpay.eq.2) then
        I_treat=2 ! treat
        I_pay=2 ! not pay
    else if (jtrpay.eq.3) then
        I_treat=1 ! not treat
        I_pay=2 ! not pay
    end if 

    

    
end if ! for being on consumption floor even if you don't pay
end if ! for having a health shock

end if

!If you don't get treated, then there is a probability the shock and R don't get reported

!initialize to not having shock
jdp_new=jdp
jdu_new=jdu
js_new=js
jr_rep=jr

!dp (2,4,6,8)
!if ((jtype2.eq.2).OR.(jtype2.eq.4).OR.(jtype2.eq.6).OR.(jtype2.eq.8) ) then
if (jdp.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,1,I_treat))   then
      jdp_new=1 
    end if 
end if 


!du (5,6,7,8)
!if ((jtype2.eq.5).OR.(jtype2.eq.6).OR.(jtype2.eq.7).OR.(jtype2.eq.8) ) then
 if (jdu.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,2,I_treat))   then
      jdu_new=1 
    end if 
end if

!s (3,4,7,8)
 if (js.eq.2) then
    rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_shocks(jh,3,I_treat))   then
      js_new=1 
    end if 
 end if
 
 ! R
 if (jr.eq.2) then
         rb=ran0(idum) ! draw a random # between 0-1
    if (rb.gt.neu_R(jo_act))   then
      jr_rep=1 
    end if 
 end if 
 
 ! but if you get a gov transfer (have Medicaid), then all shocks and R reported correctly; this should already be OK for shocks. making sure ok for R too.
 if (jtr_act.gt.0.0d0) then 
 jdp_new=jdp
 jdu_new=jdu
 js_new=js
 jr_rep=jr
 end if 
 

call TYPE2_DET(jdu_new,jdp_new,js_new,jtype2_new)

                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                        write (17,1006) jo_act !o_act_bench
                        write (20, 1005) jas_p_act !sav_bench 
                        write (18,1005)  jx_int !HC_bench 
                        write (16,1005) Z_INT_OFFER ! wage_bench ! wage offer
                        write (19,1006) jh !H
                        write (21,1006) jzt !zt_hypo
                        write (23,1006) jo
                      end if 
                      

!Write to file
 write(12,1008) jsample, jage, je, jfix, jht_orig, jmar, jh,  jr, jtype2, jtype4, jo ,   jo_act,  jzt, jx_int, jx_int_p, oop(jtype2, jh, jage,jMEcat,jo_act)*5200.0d0 , Charges(jtype2, jh, jage,jMEcat)*5200.0d0, jcons_act*5200.0d0, jassets*5200.0d0 , jas_p_act*5200.0d0,  jtr_act*5200.0d0  , wage_penalty(jo,je)*Z_INT_OFFER   , 5200.0d0*actual_earn, max(0.0d0, (hrs_offer(jo_act) -phi(je,jh,jtype2) ))*100.0d0, exp( log(wage_penalty(jo_act,je)*Z_INT_OFFER)+random_normal()*stdze(je) ) ,  inc_spouse(jage,je,jh,jtype4,jfix)*5200.0d0, jtrpay  , oop_spouse(jage,je,jtype4,jo_act)*5200.0d0 ,base_1_pay(1)*5200, base_1_pay(2)*5200, base_2_pay(1)*5200, base_2_pay(2)*5200, 0, jtype2_new, jr_rep, joption_pay, juntreated, B_sim(jtrpay)
 
 
 
!draw marital status for next period -- we do this here because it depends on jo_act up to retirement, but not after.
            rf=ran0(idum) ! draw a random # between 0-1
            
                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.40) then
            write(14,1005) rf 
            else if (jage.LE.40) then
 

        rf=ran_hypo(12+11*(jage-1))
                                end if
                                
            if (rf.lt.(Tp_mar(1,jmar,je,jage,jh,jinc,jo_act))) then
            jmar_p=1
            else 
            jmar_p=2
            end if
          jmar=jmar_p     
!Health Risk R
          rb=ran0(idum) ! draw a random # between 0-1
                                  if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.40) then
            write(14,1005) rb 
            else if (jage.LE.40) then

        rb=ran_hypo(13+11*(jage-1))
                                end if
          if (rb.lt.p_risk(jr, 1,jage,jh)) then
            jr_p=1
            else 
            jr_p=2
            end if
          jr=jr_p  
          
          
                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.6) then 
                       jr=1
                       end if
                       
              if (rb.lt.p_risk(jr_alt, 1,jage,jh_alt)) then
            jr_p_alt=1
            else 
            jr_p_alt=2
            end if
          jr_alt=jr_p_alt                    
                       
                       
!Health
          rb=ran0(idum) ! draw a random # between 0-1
                                  if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.40) then
            write(14,1005) rb 
            else if (jage.LE.40) then
 
 

        rb=ran_hypo(14+11*(jage-1))
                                end if
                                
      
                                
                                
                                
            bb=1
            prob=Tp_h(bb, jh,je,jht_orig,jage, jtype2,I_treat,jmc) 
         92 if (rb.lt.prob) then
            jh_p=bb
            else 
            bb=bb+1
                if (bb.le.n_h) then
                prob=prob+Tp_h(bb, jh,je,jht_orig,jage, jtype2,I_treat,jmc)
                go to 92
                else 
              jh_p= n_h  
                end if
            end if  
            jh=jh_p  !I'll do this after I draw the employment offer for next period - this is because i need hours in that bit of code, and hours depend on sick days which depend on health         

           !Health if no health shocks - set jtype2 to 1. 
            bb=1
            prob=Tp_h(bb, jh_alt,je,jht_orig,jage, 1,I_treat,jmc) 
        292 if (rb.lt.prob) then
            jh_alt=bb
            else 
            bb=bb+1
                if (bb.le.n_h) then
                prob=prob+Tp_h(bb, jh_alt,je,jht_orig,jage, 1,I_treat,jmc)
                go to 292
                else 
              jh_alt= n_h  
                end if
            end if 
            
            

!Next period states needed only for those of working age, but not 64: 
if (jage.NE.(n_ret-1)) then    
    
    
  
        
    if ((jExperimentNumber.EQ.0).AND.((jBenchmark_exp.EQ.23).OR.(jBenchmark_exp.EQ.28)).AND.(HC_bench(jage+1).ne.100.0d0)) then
        jx_int = HC_bench(jage+1)
    else 
           jx_int=max(0.0d0, min(jx_int_p , dble(n_hc)) )
    end if 
    


      
            
!Draw employment offer jo

      ro=ran0(idum) ! draw a random # between 0-1
                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.39) then
            write(14,1005) ro 
            else if (jage.LE.39) then

        ro=ran_hypo(15+11*(jage-1))
                                      end if
                                      
        bb=1

        prob=Pr_o_i(bb,jo_act, je,jage)

     34 if (ro.lt.prob) then
        jo=bb
        else 
        bb=bb+1
        if (bb.le.n_o) then
        prob=prob+Pr_o_i(bb,jo_act, je,jage)       
        go to 34
        else 
            jo=n_o
        end if
        end if   
        

        
        
!Draw zt next period
                               rb=ran0(idum) ! draw a random # between 0-1
                                                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.39) then
            write(14,1005) rb 
            else if (jage.LE.39) then 
        rb=ran_hypo(16+11*(jage-1))
            end if
                                      
        if (jo_act.eq.1) then
         
		            bb=1
		            prob=p_zt(je,bb,1)
89		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,1)
		            go to 89
                    else 
                        jzt=n_zt
		            end if
                    end if  
                    
 
        else 

		            bb=1
		            prob=p_zt(je,bb,2)
99		            if (rb.lt.prob) then
		            jzt=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_zt) then
		            prob=prob+p_zt(je,bb,2)
		            go to 99
                    else 
                        jzt=n_zt
		            end if
                    end if  
        end if 
        
!Spouse's work status next period (jos)- this depends on the individual's state next period, so it's ok to enter the new jh here. Also, we've already determined the new jmar next period

                        ! draw spouse work status  1=not work, 2=work ; if not married, doesn't matter, jos will be drawn but not used.
                        rb=ran0(idum) 
                                                                               if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.0.AND.jage.LE.39) then
            write(14,1005) rb 
            else if (jage.LE.39) then

        rb=ran_hypo(17+11*(jage-1))
                                      end if
                        if (rb.lt.(prob_os(jage+1,je,jh,1,jfix))) then ! this is the jh next period   
                        jos=1
                        else 
                        jos=2
                        end if 
               
end if 
end if 


		
!****************************************************************************************************************************************


if (jage.ge.n_ret) then !individual is retired	
           
!Draw health shocks jtype2
                  rf=ran0(idum) ! draw a random # between 0-1
                  
                                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                                         write(22,1005) rf 
                                       end if 
                                       
		            bb=1
		            prob=shocks_risk(je,bb, jh, jr, jage) 
		         37 if (rf.lt.prob) then
		            jtype2=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_type2) then
		            prob=prob+shocks_risk(je,bb, jh, jr, jage)
		            go to 37
                    else 
                        jtype2=n_type2
		            end if
                    end if
                    
if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.5 ) then
 jtype2=1
end if	


!Draw survival		            
                     rx=ran0(idum)
                                       if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                                         write(22,1005) rx 
                                       end if 
		             if (rx.gt.p_surv(jh,je,jmar , jage)) then !individual dies	
                         
   if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
    if (jage.eq.n_age) then
            rb=ran0(idum)
            write(22,1005) rb 
        else 
                 do jno=1, 1+(n_age-jage)*6 ! this makes sure that all individuals have 443 observations written to the random draws file- makes it easier to read.  
                     rb=ran0(idum)
            write(22,1005) rb
                 end do 
        end if 
end if                
                         
                         
                         
                     go to 33 !simulate new individual 
                     end if	             
!Draw the medical expenditure catastrophic shock
                  rf=ran0(idum) ! draw a random # between 0-1
                                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                                         write(22,1005) rf 
                                       end if 
		            bb=1
		            prob= Prob_MEcat(bb)
		         18 if (rf.lt.prob) then
		            jMEcat=bb
		            else 
		            bb=bb+1
		            if (bb.le.n_MEcat) then
		            prob=prob+Prob_MEcat(bb)
		            go to 18
                    else 
                        jMEcat=n_MEcat
		            end if
                    end if            
                     
		             
!Solve optimal consumption/sav: jcons_act, jas_p_act


!income
income = max(0.0d0, (interest-1.0d0)*jassets + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -  prem_mcare(jmar)  - max(0.d0,  oop(jtype2, jh, jage,jMEcat,1) +  oop_spouse(jage,je,jmar,1) - 0.075d0*((interest-1.0d0)*jassets + pen(jhc_ret,  je,jfix) + inc_spouse(jage,je,1,jmar,jfix) ) )) !taxable income. the last term is OOP in excess of a threshold of your income, which is deductable.
base_2 = prem_mcare(jmar) + oop(jtype2, jh, jage,jMEcat,1) + oop_spouse(jage,je,jmar,1)  + max(0.0d0,  (( (income*5200.0d0) - tax_lambda(jmar) *( (income*5200.0d0)**(1.0d0-tax_tau(jmar))) )/5200 )) ! total resources subtracted: premiums, OOP, and tax, at the HH level

                
        
!Initialize
 RET_act=0.0d0
 jtr_act=0.0d0
 jcons_act=0.0d0
 jas_p_act=0.0d0
 maxB=0.0d0
  
if (((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) .le. (cons_floor_at(jage, jmar,je,jh)) ) then    ! get transfer                                                                                                                                  
    jtr_act= cons_floor_at(jage, jmar,je,jh) - (interest*jassets+ pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2) 
    jcons_act=cons_floor(jage, jmar,je,jh)
    jas_p_act=0.0d0

else
   jtr_act=0.d0
!Find max savings point
   MAX_SAV=  interest*jassets  + pen(jhc_ret, je,jfix) + inc_spouse(jage,je,1,jmar,jfix) -base_2  -cons_floor_at(jage, jmar,je,jh) -jassets
   
!Find lower bound
            jsav_upper_ret=0
          do while ( MAX_SAV .ge. sav_sim(jsav_upper_ret+1) )
           jsav_upper_ret=jsav_upper_ret +1
           if (jsav_upper_ret.EQ.n_sav_sim) exit
          end do 
          
               if (jsav_upper_ret.lt.1) then
                 PRINT *,'Error'
               end if 
         
!Find value of consuming all and saving nothing
   CON=((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) /(1.0d0+taxc)      
        call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
                        V_death =  Bequest_val(0.0d0)  + cost_death
         !Calculate expected value associated with MAX_SAV assets tomorrow
         
                if (jage.eq.n_age) then
                 V_c_all_trpay=UTI+ beta_hat(je)* V_death  
                else 
                   EXPECTED_RET_PRIME=0.0d0
                    do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1, n_MEcat; do jmar_p=1, n_mar
                     EXPECTED_B_PRIME = EXPECTED_B_PRIME + Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p)* Tp_mar(jmar_p,jmar,je,jage,jh,jinc,jo_act) * (RET_PRIME(jhc_ret,jMEcat_p,1,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p) * p_surv(jh_p,je,jmar_p, jage+1) +  V_death* (1.0-p_surv(jh_p,je,jmar_p, jage+1)) )                                      
                    end do; end do; end do  ; end do; end do 
                 
                 V_c_all_trpay = UTI + beta_hat(je)* EXPECTED_B_PRIME 
                end if 
         
!Find value of consuming the consumption floor and saving as much as possible 
      call func_uti_ret(je,jh,jmar,jtype2,cons_floor(jage, jmar,je,jh),UTI)
      ASP= MAX_SAV+jassets
      V_death =  Bequest_val(MAX_SAV+jassets)  + cost_death    
                   
                if (jage.eq.n_age) then
                 V_cf_trpay=UTI+ beta_hat(je)* V_death  
                else 
                  EXPECTED_RET_PRIME=0.0d0   
                do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1,n_MEcat; do jmar_p=1,n_mar
                    call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) *  Tp_mar(jmar_p,jmar,je,jage,1,1,1) * ( p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) + V_death* (1.0 -p_surv(jh_p,je,jmar_p, jage+1) ) )
                end do; end do; end do  ; end do; end do
        V_cf_trpay=UTI+ beta_hat(je)*EXPECTED_RET_PRIME
                end if 
                 
               

!loop for optimal savings
         maxB =V_cf_trpay
         jas_p_act=MAX_SAV+jassets
         jcons_act=cons_floor(jage, jmar,je,jh)
   do jsav= jsav_upper_ret , 1, -1
       
!Find CON, ASP and UTI

    ASP=  jassets+  sav_sim(jsav)
                if (ASP.lt.0.0d0)  exit
                V_death=cost_death + Bequest_val(ASP)
                
        CON= ( ((interest-1.0d0))*jassets  + pen(jhc_ret,  je,jfix)+ inc_spouse(jage,je,1,jmar,jfix) -base_2 -sav_sim(jsav) )/(1.0d0+taxc)

        if (CON.lt.cons_floor(jage, jmar,je,jh)) then
         PRINT *,'Error'
          pause
        end if      
 
     call func_uti_ret(je,jh,jmar,jtype2,CON,UTI)
     
!Find continuation value
            
                EXPECTED_RET_PRIME=0.0d0 
                
                if (jage.eq.n_age) then
                 RET_act=UTI+ beta_hat(je)* V_death  
                else 
                    
                do jtype2_p=1, n_type2 ; do jh_p=1, n_h ; do jr_p=1, n_r ; do jMEcat_p=1,n_MEcat; do jmar_p=1,n_mar
                    call InterpolateScalar(as(1:n_as,jtype1,jage+1), RET_PRIME(jhc_ret,jMEcat_p,1:n_as,jtype1,jh_p,jr_p,jtype2_p,jage+1,jmar_p), ASP, RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) )
                EXPECTED_RET_PRIME= EXPECTED_RET_PRIME+ Tp_h(jh_p, jh,je,jht_new,jage, jtype2,2,1)*p_risk(jr, jr_p,jage,jh)*shocks_risk(je,jtype2_p,jh_p,jr_p, jage+1)*Prob_MEcat(jMEcat_p) *  Tp_mar(jmar_p,jmar,je,jage,1,1,1) * ( p_surv(jh_p,je,jmar_p, jage+1)* RET_PRIME_INT(jMEcat_p,jh_p,jr_p,jtype2_p,jmar_p) + V_death* (1.0 -p_surv(jh_p,je,jmar_p, jage+1) ) )
                end do; end do; end do  ; end do; end do
        RET_act=UTI+ beta_hat(je)*EXPECTED_RET_PRIME
                end if 
                
!Find max value
                if (RET_act.lt.maxB) exit
                maxB=RET_act 
                jas_p_act=ASP
                jcons_act=CON         
    		    	    
56 end do !for jsav for retired ppl
   
     if  (V_c_all_trpay.gt.maxB) then
      maxB =V_c_all_trpay
      jas_p_act=0.0d0
      jcons_act= ((interest*jassets+ pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 )) /(1.0d0+taxc)
     end if 
     
end if ! for being on tr
      


   
!Write decisions to file
write(12,1008)  jsample, jage, je, jfix, jht_orig, jmar, jh,  jr, jtype2, 0, 0 ,  0,  0, 0.0d0 , 0.0d0, oop(jtype2, jh, jage,jMEcat,1)*5200.0d0 , 0.0d0 , jcons_act*5200.0d0, jassets*5200.0d0 , jas_p_act*5200.0d0,  jtr_act*5200.0d0  , 0.0d0   , pen(jhc_ret,  je,jfix)*5200.0d0 ,  0.0d0 ,  0.0d0,  inc_spouse(jage,je,1,jmar,jfix)*5200.0d0, 0  , oop_spouse(jage,je,jmar,1)*5200.0d0, ((interest*jassets+  pen(jhc_ret,  je,jfix)  + inc_spouse(jage,je,1,jmar,jfix) - base_2 ))*5200.0d0 , 0.0, base_2*5200.0d0, 0.0, jhc_ret, 0, 0, 0, 0, maxB

if (jage.eq.n_age) then
go to 33 !simulate new individual      
end if                            
      
! draw marital status for next period
            rf=ran0(idum) ! draw a random # between 0-1
                                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                                         write(22,1005) rf 
                                       end if 
            if (rf.lt.(Tp_mar(1,jmar,je,jage,1,1,1))) then
            jmar_p=1
            else 
            jmar_p=2
            end if
          jmar=jmar_p            
!Risk factors
          rb=ran0(idum) ! draw a random # between 0-1
                                        if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                                         write(22,1005) rb 
                                       end if 
            bb=1
            prob=p_risk(jr, bb,jage,jh)
         39 if (rb.lt.prob) then
            jr_p=bb
            else 
            bb=bb+1
            if (bb.LE.n_r) then
            prob=prob+p_risk(jr, bb,jage,jh)
            go to 39
            else 
            jr_p=n_r  
            end if
            end if
           jr=jr_p        
!Health
          rb=ran0(idum) ! draw a random # between 0-1
                                      if (jExperimentNumber.EQ.0.AND.jBenchmark_exp.EQ.21) then
                                         write(22,1005) rb 
                                       end if 
            bb=1
            prob=Tp_h(bb, jh,je,jht_orig,jage, jtype2,2,1) ! YOU ALWAYS GET TREATED IF RETIRED
         32 if (rb.lt.prob) then
            jh_p=bb
            else 
            bb=bb+1
            if (bb.le.n_h) then
            prob=prob+Tp_h(bb, jh,je,jht_orig,jage, jtype2,2,1)
            go to 32
            else 
            jh_p= n_h 
            end if
            end if 
    	    jh=jh_p          
end if !for being retired
           
!****************************************************************************************************************************************            
!Assets next period  
  if (jage.lt.n_age)then  
            jassets=jas_p_act          
  end if                 		            
!****************************************************************************************************************************************
  
end do !age
33 end do !sample
   end do; end do ! jht, jfix

close (12)
close (13)
close (14)
close (15)
close (16)
close (17)
close (18)
close (19)
close (20)
close (21)
close (22)
close (23)
end do
End Subroutine